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Clumpy Outflow from Supercritical Accretion Flows 
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Abstract 

Significant fraction of matter in supercritical (or super-Eddington) accretion flow is blown away by radiation 
force, thus forming outflows, however, the properties of such radiation-driven outflows have been poorly under- 
stood. We have performed global two-dimensional radiaion-magnetohydrodynamic simulations of supercritical ac- 
cretion flow onto a black hole with 10 or IO^Mq in a large simulation box of 514rs x 514rg (with rs being the 
Schwarzschild radius). We confirm that uncoUimated outflows with velocities of 10 percents of the speed of light 
emerge from the innermost part of the accretion flow over wide angles of 10°-50° from the disk rotation axis. 
Importantly, the outflows exhibit clumpy structure above heights of ^ 250rg. The typical size of the clumps is 
^ lOrs, which corresponds to one optical depth, and their shapes are slightly elongated along the outflow direction. 
Since clumps start to form in the layer above which (upward) radiation force overcomes (downward) gravity force, 
Rayleigh-Taylor instability seems to be of primary cause. In addition, a radiation hydrodynamic instability, which 
arises when radiation funnels through radiation-pressure supported atmosphere, may also help forming clumps of 
one optical depth. Magnetic photon bubble instability seems not to be essential, since similar clumpy outflow struc- 
ture is obtained in non-magnetic radiation-hydrodynamic simulations. Since the spatial covering factor of the clumps 
is estimated to be ~ 0.3 and since they are marginally optically thick, they will explain at least some of rapid light 
variations of active galactic nuclei. We further discuss a possibility of producing broad-line clouds by the clumpy 
outflow. 

Key words: accretion, accretion disks — instabiUties — radiative transfer — ISM: clouds — ISM: jets and 
outflows 



1. Introduction 

There is growing observational evidence supporting that sig- 
nificant amount of material is blown away from accretion flow 
(or disk) onto black holes in forms of outflows and jets. In most 
cases, the outflows are observed through blue-shifted absorp- 
tion line features (e.g., MiUer et al. 2004; Kotani et al. 2006; 
Kubota et al. 2007; Neilsen et al. 2012; Ponti et al. 2012 for 
Galactic sources and Terashima & Wilson 2001; Pounds et al. 
2003; Reeves et al. 2003; Ganguly & Brotherton 2008; Pounds 
& Reeves 2009 for active galactic nuclei: AGNs). The esti- 
mated mass outflow rate sometimes exceeds the mass accretion 
rate to the central black holes (e.g.. Miller et al. 2008; Ueda et 
al. 2009; Neilsen et al. 2011). The dense outflows can also 
produce Compton upscattering soft photons from the underly- 
ing accretion flow, thereby making the emergent spectra harder 
(e.g, Gladstone et al. 2009; Kawashima et al. 2012). The struc- 
ture of the accretion flow is also affected by the presence of 
outflows (Kawabata & Mineshige 2009). 

In addition, the outflows are expected to produce deep im- 
pacts on their environments, since the outflows from black hole 
accretion flow, especially those from the innermost part of the 
flow, do carry large amount of mass, momentum, and energy. 
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Such powerful outflows may cause dynamical feedback which 
may trigger (or quench) star formation in surrounding inter- 
stellar and galactic medium, thus promoting metal enrichment 
(Murray et al. 1995; Granato et al. 2004; Di Matteo et al. 2005). 
Larger impacts are expected by higher luminosity objects with 
larger mass accretion rates. 

At luminosities close to the Eddington luminosity, radiation 
energy density overcomes matter energy density so that radia- 
tion field can also be dynamically important; radiation-driven 
outflow thus naturally arises. Matter emits, absorbs, and scat- 
ters radiation, while radiation gives (or removes) energy and 
momentum to matter In this way, matter and radiation are 
strongly coupled to each other. To study the outflow from lu- 
minous systems, we thus need to properly simulate strong cou- 
pling between matter and radiation. 

There are two basic lines of simulation study of outflow from 
luminous black hole accretion flow: hydrodynamic simulations 
of line-driven outflow by Proga and collaborators (Proga et al. 
1998; Proga et al. 2000), and radiation-hydrodynamic simu- 
lations of continuum-driven outflow by our group (Ohsuga et 
al. 2005; Takeuchi et al. 2009). By performing the hydrody- 
namic simulations incorporating radiation force due to spec- 
tral lines (line force) into the equation of motion, Proga & 
Kallman (2004) have studied the dynamics of the line-driven 
outflow. They investigated the outflow properties above the 
disk by setting an optically thick and geometrically thin disk 
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on the equatorial plane as a boundary condition. The radiation 
energy equation is not solved and the reprocessed photons are 
neglected. In contrast, our simulations self-consistently solve 
the accretion disks and the outflows, although our studies are 
restricted to the continuum-driven outflows. By solving radia- 
tion energy equation, we consider diffusive photons which suf- 
fer numerous Thomson scattering. The continuum radiation 
force via diffusive photons drives outflows. 

In traditional hydrodynamic (or radiation-hydrodynamic) 
models of accretion flows, the disk viscosity, most important 
physical process in the disk theory, is described by the phe- 
nomenological a-viscosity model, although it should, in prin- 
ciple, be determined in accordance with magnetic field ampli- 
fications and dissipation within the flow. In other words, global 
multi-dimensional radiaion-magnetohydrodynamic (radiation- 
MHD) simulations are necessary for understanding luminous 
black hole inflow-outflow system. Under such considerations, 
we have recently started global two-dimensional radiation- 
MHD simulations (Ohsuga et al. 2009; Ohsuga & Mineshige 
2011, hereafter Papers I and II, respectively). The most im- 
portant findings in those simulations are ubiquitous outflow in 
every mode of accretion flow and a new type of jet emerging at 
high luminosities; magnetically collimated, radiation-pressure 
driven jet (Takeuchi et al. 2010). 

The aim of the present study is to examine large-scale be- 
haviour of continuum-driven outflow from supercritical (or 
super-Eddington) accretion flows. For this purpose, we run 
new simulation with a significantly increased simulation box 
size of 514rs x 514rs, where rs is the Schwarzschild radius 
(Note that it was only 105rs x 103rs in Papers II). With this 
treatment we could find clumpy outflow structure in distant re- 
gions from the central black hole. The plan of this paper is as 
follows. In the next section, we overview the basic equations 
of our radiation-MHD simulations. We present the simulation 
results of clumpy outflow from supercritical accretion flows in 
section 3. The formation mechanisms of clumpy structure and 
the observational implications will be also discussed in section 
4. The final section is devoted to concluding remarks. 

2. Basic Equations 

We simulated supercritical accretion flows and the associ- 
ated outflows by using the global two-dimensional radiation- 
MHD code developed by Papers I and II. The basic equations 
and the numerical method are described them in details. Hence, 
we overview them and stress the differences from the previous 
calculations. 

We used cylindrical coordinates {R, 9, z), where R, 9, and 
z are the radial distance, the azimuthal angle, and the verti- 
cal distance, respectively. We assume that the flow is non- 
self-gravitating, reflection symmetric relative to the equatorial 
plane (with z = 0), and axisymmetric with respect to the ro- 
tation axis (i.e., d/d9 = 0). The basic equations which take 
the terms up to the order of (v/c)^ are as follows (Mihalas & 
Weibel Mihalas 1984; Stone & Norman 1992): the continuity 
equation, 

dp 
dt 



V • (pv) - 0, 



(1) 



the equations of motion, 

d{pv) ^ f BB 
^ ' ' V • pvv - 



dt 



- -Fo -joVi/'PN, 
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the energy equation of the gas, 

de 



dt 



-V-(eu) = -pgasV-t>- 



47r 



- 47rXabs-Bbb + CXabs-Bo, 

the energy equation of the radiation, 

dEa 



dt 



- V • {Eqv) = - V • Fo - Vt; : Po 

+ 47rXabsBbb - CXabs^^O, 



and the induction equation, 

dB ^ f „ 47r ^ 
— = Vx lvxB-—rjJ 



(3) 



(4) 



(5) 



Here, p is the matter density, v is the flow velocity, c is the 
speed of light, e is the internal energy density of the gas, pgas 
is the gas pressure, B is the magnetic field, J = c V x S/47r is 
the electric current, Bbb = f^T^ns/^ the blackbody intensity, 
(T is the Stefan-Boltzmann constant, Tgas is the temperature of 
the gas, Eq is the radiation energy density, Fq is the radiative 
flux, Pq is the radiation-pressure tensor, Xabs = {^s + Khi)p 
is the absorption opacity per unit volume (with dimensions of 
length"^), X = ('^es + Kff + Kbf )p is the total opacity per unit 
volume, and rj is the resistivity, respectively. The subscript 
means the values measured in the co-moving (fluid) frame. For 
simplicity, we adopt the gray (frequency-integrated) approxi- 
mation for the radiation terms. 

We adopt the pseudo-Newtonian potential, i/'pn, to incorpo- 
rate the general relativistic effects, given by -0pn = —GM/ (r — 
rs) (Paczynski & Wiita 1980). Here, r = (i?^ + z^^^ is the 
distance from the origin. The Schwarzschild radius is given by 
rs = 2GM/c^, where G is the gravitational constant, and M is 
the mass of the black hole, respectively. 

The set of equations (l)-(5) can be closed by using an ideal 
gas equation of state, pgas = (7 — l)e = pfceTgas/Af^P' 
by adopting the flux-limited diffusion (FLD) approximation to 
evaluate Fq and Pq (Levermore & Pomraning 1981). Here, 7 
is the specific heat ratio, fee is the Boltzmann constant, p is the 
mean molecular weight, and nip is the proton mass. 

We consider the electron scattering, Kcs, the Rosseland mean 
free-free absorption, ks, and bound-free absorption opacity, 
Kbf , for solar metallicity; 
-1 



Kff = 1.7x10 



-25^-2 



-7/2 



and 



: 4.8 X IQ-^'^m^^pT 
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where <tt is the Thomson scattering cross-section. 

We initially set a magnetized rotating torus surrounding a 
non-rotating black hole at the origin. The center of the torus, 
where the matter density is at maximum, is located at i? = 
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40 rg. The initial, maximum density (at the center of the torus) 
is the same as that of Model A of Paper I and II; i.e., po = 1 g 
cm~'^. We assume initially closed poloidal magnetic field and 
their strengths are determined so as to achieve the plasma-/? 
(=Pgas/Pmag) to be 100 withiu the torus, where Pmag = S^/Stt 
is the magnetic pressure. This torus is sandwiched by a non- 
rotating, non-magnetized isothermal corona. Since the initial 
ambient gas is finally ejected out of the computational domain, 
it does not affect the resulting flow structure. 

We adopt free boundary conditions at the upper and the outer 
boundaries, where we allow mass, momentum, and energy to 
free go out but nothing can enter the calculation box. At the 
inner boundary at i? = 2rs around the black hole we remove 
all the mass, momentum, and energy which enter the inner zone 
inside this radius. Technically, we first solve the non-radiative 
MHD equations for initial 1 s, and then turn on the radiation 
terms. 

The calculation methods, the initial model, and the bound- 
ary conditions so far described are the same as those adopted in 
Papers I and II. Only a difference is the size of the calculation 
box. Since the purpose of this study is to examine global gas 
dynamics of the outflow originating from supercritical accre- 
tion flow, we need a wider spatial range. Thus, we set the com- 
putational domain of cylindrical shells of 2rs < i? < 514rs and 
< z < 514rs. This large calculation box leads us to a discover 
of clumpy outflow structure, as we see in the next section. The 
grid spacing of the radial distance, AR, and the vertical dis- 
tance, Az, are uniform, given by AR = Az = 0.4rs. Note that 
the computational domain of Papers II was 2rs < R< 105 rs, 
and < z < 103rs, and the grid spacing was AR ~ Az = 0.2rs. 

3. Properties of Clumpy Outflow 

3.1. Overview of Clumpy Outflow 

Let us first overview the global gas dynamics. Figure 1 
shows the two-dimensional distributions of the matter density 
(upper panels) and the ratio of the gas temperature to the ra- 
diation temperature (lower panels) in the whole computational 
domain in the left panels and in a clumpy outflow region (de- 
fined later) in the right panels at the elapsed time of f = 9 s. The 
black hole mass is set to be M = IOA/q. The radiation temper- 
ature is given by Trad = (cEo/iay^'^. Yellow arrows in the up- 
per left panel indicate the velocity field whose values are larger 
than the escape velocity, Vcsc = (2GA//r)^/^. Green lines in 
the upper left panel indicate the surface, where the equality, 
X-Fo/c = pIV-^pnI, holds in time average between t = 8.8-9.2 
s. Upward radiation force overcomes downward gravity force 
of the central black hole in the upper and right region. 

We evaluated the following basic quantities: the mass ac- 
cretion rate is Mace ~ IOOLe/c'^; the mass outflow rate is 
-^out ~ IOLe/c'^; the photon luminosity is L ^ Le- Here, the 
Eddington luminosity Le is given by Le = AncGM / Kcs- The 
mass accretion rate and the mass outflow rate were calculated 
by summing up the mass passing through the inner boundary 
and the upper boundary per unit time with higher velocities 
than the escape velocity. The photon luminosity were calcu- 
lated based on the radiative flux at the upper boundary: 



2rs 



Mace = -2 / 2'KRpVRdz, 

Jo 

Mout = 2 / 2nRpv,dR, 

J Irs 



and 



(•514rs 

L = 2 / 2t: R{F^ + v,_EQ)dR. 



(9) 
(10) 

(11) 



The mass accretion rate onto the black hole exceeds the critical 
rate giving rise to the Eddington luminosity; i.e., the simulated 
flow is supercritical (Shakura & Sunyaev 1973; see chap. 10 of 
Kato et al. 2008 for a review). The inflow-ouflow structure in 
the central region (i?^ lOOrs, z ^ 1007's) is consistent with the 
previous result (Paper II). Note that we adjusted the color scale 
in such a way to clearly show the flow structure in the distant 
region from the central black hole. As a result, the color scale 
in the central accretion flow region is saturated. 

Figure 1 shows clumpy structure in the distant outflow re- 
gion {R ^ 2007's and z > 250rs), in which upward radiation 
force overcomes downward gravity force. Gas clumps whose 
matter density is about 10^^ g cm^'^ and velocity is about 10% 
of the speed of light are blown away over wide angles of 10°- 
50° from the disk rotation axis. 

Besides the uncollimated outflows we also confirmed the 
ejection of a high-speed, collimated jet in a narrow region 
along the disk rotation axis. This high-speed outflow is acceler- 
ated by continuum radiation force while it is collimated by the 
Lorentz force of a magnetic tower structure despite radiation 
energy greatly dominating the magnetic energy (see Takeuchi 
et al. 2010). Note that the magnetic tower structure is created 
by the inflation of toroidal magnetic field which are accumu- 
lating around the black hole. In the present paper, we focus 
properties of the clumpy outflow. 

We show the clumpy structure in more details in the upper 
right panel of Figure 1 . The depicted area is 3007's <R< 400rs 
and 350rg < z < 450rs (hereafter called as the "clumpy out- 
flow region"), corresponding to the location of the white square 
in the upper left panel. Clearly displayed is grossly inhomo- 
geneous density pattern which is composed of high density 
clumps (pci ^ 10^^ g cm^'^) and ambient low-density media 
(Pamb ^ 10~* g cm-3). 

The lower right panel shows the ratio of the gas temperature 
to the radiation temperature in the clumpy outflow region. It is 
interesting to note that the radiative equilibrium (Qjtj^j ^ Q^d' 
and hence Tgas ^ T^ad) is achieved within clumps, as well as 
within the underlying accretion flow. Here, Qjt^d ^ cxabs^^o 
and Qrad ^ 47rXabsSbb are the radiation heating rate and the ra- 
diation cooling rate, respectively. In the ambient media, in con- 
trast, decoupling of matter and radiation occurs (Tgas > Tlad) 
because of low matter density. 

Also note that Figure 1 gives a snapshot. In fact, the clumpy 
structure is not stationary but changes its shape in time. The 
gas particles that compose a clump are not the same. As we 
will see later the gas particles move faster than clumps and 
sometimes go through a clump. 
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Fig. 1. Two-dimensional stmcture of a supercritical accretion flow and the associated outflow around a black hole with M = 10A/q . The upper and 
lower left panels, respectively, show the color contours of the matter density overlaid with the flow velocity vector and those of the ratio of the gas 
temperature to the radiation temperature in the whole computational domain. The velocity arrows are displayed only in the region in which their values 
exceed the escape velocity. Green Unes in the upper left panel indicate the surface which upward radiation force equals downward gravity force of the 
central black hole. The right two panels are the same as those in the left panels but in a narrower region; SOOrs < < 400rs and 350rs <2:<450rs. 
The elapsed time is 9 s, which corresponds to about 30 times Keplerian timescale al R = AOr<j . 



3.2. Clump size and shape 

The two-dimensional contour plots shown in the previous 
section are useful to understand the clumpy patterns of the 
outflow, but it is not always easy to derive a typical size of 
the clumps and a mean separation between the clumps. For 
deriving such statistical quantities, auto-correlation function 
(ACF) is quite useful (see Appendix for the actual calculation 
method). 

We show ACFs of the matter density as a function of the typ- 
ical height at the elapsed time of t = 9 s in Figure 2. The inte- 
gration range of the ACFs are (i?in,i?out) = (100rs,200rs) for 
z = 160rs, (i?i„,i?out) = (150rs,250rs) for z = 230-250rg, 
and (i?in,^out) = (300rg,400rg) for z — SSOrg, respectively. 

The ACFs show several noteworthy features. First, no strong 
correlation is found at low altitudes, z ^ 200rs, meaning that 
no clumps there are. Second, we see a sharp peak in the ACFs 



for the data above z ^ 250rg, indicating that clumps are being 
formed at ^ 250rg. The typical clump size tc\ can be measured 
by the width of the central peak and is ^ lOrg. Third, the shape 
of ACFs at higher z ^ 250rs do not change appreciable. This 
means, the clumps retain their shapes in the statistical sense 
even after moving upwards. The clump size is 25 times the 
grid spacing (Ar — Az = 0.4rg). We can thus conclude that 
each clump is well resolved in our simulations. 

To understand what the typical size of ~ lOrg means phys- 
ically, we estimate the optical depth of clumps for this length 
scale inserting typical clump density, pc\ ^ 10~^ g cm~^, find- 
ing 

Tcl = KesPcl4l ^ 1- (12) 

That is, the clump size is on the order of one optical depth. This 
is a very important property of the clumps to consider their for- 
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Fig. 2. Auto-correlation functions (ACFs) of the matter density in the 
/J-direction (<5_R) as a function of the typical height at the elapsed time 
of 4 = 9 s. The width of the central peak represents the typical clump 
size, whereas the separation to its neighbouring peak represents the 
typical clump separation. 

mation mechanism. Here, we note that the electron scattering 
opacity dominates over the absorption (free-free and bound- 
free) opacities because of high gas temperature in the present 
simulation; that is, xl P ^ ^cs- This relation quite generally 
holds in radiation-dominated accretion disks or flows. 

We here integrated ACFs in the i?-direction. Note that 
clumps constitute an elongate shape along the outflow direc- 
tion, which is obvious in the upper panels of Figure 1 . This 
feature is also important to understand their physical formation 
mechanism. 

5.5. Anti-correlation between matter density and radiation 
force 

To understand the physical cause of the clump formation, we 
need to know the relationship between different physical quan- 
tities. The fact that the typical clump size corresponds to one 
optical depth indicates the relevant physical mechanism under- 
lying the clump formation to be somehow related to radiation 
processes. 

Let us check the magnitude and the direction of the radiation 
force per unit mass, /rad = X-Fo / pc, in the contour plot of the 
matter density, which is shown in Figure 3. Anti-correlation 
between the matter density and the absolute magnitudes of the 
radiation force is clear in the clumpy outflow region; i.e., we 
see longer arrows (representing stronger radiation force) in the 
low-density regions (with dark colors). It looks as if the radia- 
tive flux avoids dense (clumpy) regions and instead selectively 
pass through low-density channels between the clumps. The 
direction of the radiative flux seems to be responsible for the 
elongate shape of clumps. 

To see this anti-correlation nature more explicitly, we plot 
in Figure 4 the values of matter density and radiation force in 
each grid point in the clumpy outflow region. Again, the anti- 
correlation is clear Moreover, the relationships between theses 
two quantities are roughly expressed by a power-law, 

/radCXp-^ (13) 



log p [g/cm^] 




R/r, 

Fig. 3. Color contour of matter density overlaid with radiation force 
vector per unit mass (arrows) in the clumpy outflow region (correspond- 
ing to the region depicted in the right panels of Figure 1). The arrows 
are displayed only in the region where their absolute values exceed 
5 X 10'-' dyng-i. 

where a power-law index is found to be at 10~^ 

g cm^'^ while a ~' at p ^ 10^^ g cm^'^. The break of the 
power-law index is understood in terms of the optical depth. In 
the optically thick media, where the radiative diffusion approx- 
imation holds, the radiation flux is inversely proportional to the 
optical depth (and, hence, the matter density); i.e., 

/radOcFoOC^. (14) 

p 

In the optically thin media where the streaming limit applies, 
on the other hand, the absolute value of the radiation force is 
given by 

/rad oc Fo OC E'o. (15) 

We thus conclude that the break point of the power-law in- 
dex indicates the boundary between optically thick and thin 
limit. As long as the radiation energy density i^o is more 
or less homogeneous, which is actually the case, the radia- 
tion force is weaker in the optically thick medium, roughly in- 
versely proportional to the optical depth, /lad oc r^^. In the 
perturbed radiation-supported atmosphere, the radiative flux 
tends to avoid high-density regions. 

Let us illustrate the anti-correlation in a somewhat different 
way. Figure 5 shows the cross-correlation functions (CCFs) 
between the matter density and the absolute magnitudes of ra- 
diation force per unit mass. The integration range is the same 
as Figure 2. Anti-correlation is expressed by the negative value 
at zero separation. We also notice that these anti-correlations 
actually evolve as the flow moves upward (towards higher z). 

3.4. Structural Variations of a Gas Element 

In order to understand how clumps form and grow at z ^ 
250rs, we pick up one test gas element, follow its trajectory, 
and plot the time evolution of its physical quantities. The se- 
lected gas elements travel from {R, z) = (94rs, ISOrs) (within 
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Various timescales Definitions Values 



Dynamical time scale of outflow 


200rs/t; 


> 


10- 


-^s 


Sound crossing time scale 


4l/Cs=4lbgas/p)-'/' 


> 


10- 


-Is 


Joule heating time scale 




> 


10- 


-2 s 


Radiative heating time scale 


e/Qid = e/cXabs-Bo 


> 


10- 


-4 s 


Radiative cooling time scale 


e/Qrad = e/47rXabs-Bbb 


> 


10- 
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Fig. 4. Correlation diagram between matter density and the absolute 
magnitudes of radiation force per unit mass in the clumpy outflow re- 
gion (corresponding to the area of the right panels of Figure 1). 

the accretion flow) at the elapsed time of t = 8.8 s and reaches 
the point {R,z) = (390rs, 480rs) at i = 9.7 s. The trajectory 
of the element is calculated by the following equation, 

r"-^i =r" + 7;"Af, (16) 

where r" (n = 1 , 2, • • •) is the position vector of the test particle 
at the n-th time step (t = <"), is the velocity vector att~ t" 
obtained from the simulation data, and At — t"-''^ — is the 
time interval. We set the time interval to be At ~ 10^"', which 
is shorter than any physical time scales (Table 1). 

Figure 6 show a typical time evolution of gas elements. 
Here, the top, middle, and bottom panels indicate gas and ra- 
diation temperatures, the ratio of the radiation energy density 
to the sum of the gas and magnetic energy densities expressed 
by £ = Eq / {e + B'^ /Stt), and the matter density, respectively. 
We find that most gas elements go through the low gas density 
and the high gas temperature region before forming the clump. 
Initially (at t ^ 9.2 s), this gas element was within the dense 
accretion flow, in which matter and radiation is strongly cou- 
pled; i.e., Tgas = Ilad. Alt ^ 9.2-9.5 s, the decoupling occurs 
as the matter density decreases. Finally, the hot gas element 
merges with pre-existed cold clump and cools down (p 10~® 
g cm-3; Tgas - 10'' K). 

The gas temperature in the hot region (inter-clump space) 
is much larger than the radiation temperature due probably to 
compressional heating caused by weak shock. We have con- 
firmed that kinematic energy density is comparable or slightly 
exceeds the thermal energy of gas. The coupling between the 
radiation and matter is so weak there that the radiative cooling 



Fig. 5. Cross-con'elation functions (CCFs) between the matter density 
and the absolute magnitudes of radiation force per unit mass in the 
i?-direction (SR) as a function of the typical height at the elapsed time 
of t = 9 s. The negative value at the zero separation means the anti- 
correlation between each quantity. 

cannot be effective. (This contrasts with the case of clump re- 
gion, in which tight coupling between the radiation and matter 
leads to the equal temperatures.) Further, Joule heating can 
assist the energy gain but cannot be a main heating source, 
since we will show in section 4.2 that similar high tempera- 
tures are obtained by radiation-MHD simulations with no mag- 
netic fields. The gas entropy is much higher in the hot region 
than that within the clumps. Therefore, the hot region is buoy- 
ant. Cautions should be taken here, however, since inverse- 
Compton scattering is not considered in the present simula- 
tions. We expect somewhat lower temperature in the inter-blob 
space (Kawashima et al. 2009). 

To see what force is responsible for the acceleration of the 
clumpy outflow, we compare the time-averaged forces over t = 
9.60-9.65 s by using the data in Figure 6 (see Table 2). Here, 
/con = v^R^^ is the centrifugal force, /Lor = -p"^ Vpmag + 
p~^V • {BB /An) is the Lorenz force (the magnetic -pressure 
force plus the tension force), /gas = — p~^Vpgas is the gas- 
pressure force, cos6c\ = vf!{v\ + v1)~^^'^ wAsm9c\ = Vz{v\ + 
w^)~^/2 are the radial and vertical components of the acceler- 
ation of the gas particle, respectively. We understand that the 
continuum radiation force is the driving force of the clumpy 
outflow. This result is reasonable, since the clumpy outflow 
region is totally radiation-energy dominated (see the middle 
panel in Figure 6). 



No. 



Clumpy Outflow Formation 

Table 2. Comparison of the forces asserted on the clumpy outflow. 
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Fig. 6. Time variations of some physical quantities of one test clump 
element. Top panel: the gas temperature (solid lines) and the radia- 
tion temperature (dashed line). Middle panel: the ratio of the radiation 
energy density to the sum of the gas and magnetic energy densities. 
Bottom panel: the matter density. This gas element emerging from the 
accretion flow at the elapsed time of t = 8.8 s, first moves through a 
low density, high temperature inter-clump region region at t ~ 9.2-9.5 
s, and merges into a high density, low temperature clump at t ~ 9.6 s. 

4. Discussion 

4.1. Mechanism of creating clumpy outflow 

As we have seen in the previous section, the clumpy outflow 
has a number of unique features: (1) clumpy structure appears 
in the layer where upward radiation force overcomes down- 
ward gravity force, (2) a clump size is about one optical depth, 
(3) there is a clear anti-correlation between the matter density 
and radiation force, (4) temperature variations of some clumps 
are not monotonic increase nor monotonic decrease. On the 
basis of these facts we discuss plausible physical mechanism 
of clump outflow. 

The fact (4) obviously indicates that clump formation cannot 
be explained by a thermal instability, which cause that the gas 
temperature monotonically decreases (increases) with a mono- 
tonic increase (decrease) of the matter density (e.g., Balbus & 
Soker 1989). The cooling rate is expressed as Q~^^ oc p^Tgls oc 

3/2 

Tgas under the condition of constant gas-pressure, which is a 
good assumption in the present simulations (see Fig. 6). Since 
the thermal instability criterion is written as (e.g. Kato et al. 
2008) dQ^^^/ dTgas > dQ^^^/ dTgas, the system should be un- 
stable, if the heating rate ^'^^^ not depend on temperature. 
This is the case corresponding to the interstellar medium. In 
the present case, we evaluate oc (sff + km)pEo Tgas^^^ 



— 7/2 

for the Kramers-type opacity (ng + /Cbf oc pTgas ) adopted 
in the present study and under the conditions of constant gas- 
pressure and constant radiation-energy density (Eq). Then, the 
system should be thermally stable. However, it is not easy to 
express the heating rate in such simple forms for such dynam- 
ically evolving media as those we study, since compressional 
heating may play a key role. 

The fact (1) means that the clumpy outflow is Rayleigh- 
Taylor unstable. In luminous black hole inflow-outflow sys- 
tem, radiation field acts like an effective gravitational field. 
The green line of the upper left panel in Figure 1 indicates that 
the direction of the effective gravitational field is reversed by 
the strong radiation force of the flow. Since the matter density 
decreases in the direction of the acceleration, Rayleigh-Taylor 
instability should set out. The overturning motion of gas ele- 
ments will form seeds of clumpy density pattern. We thus con- 
clude that Rayleigh-Taylor instability is the primary cause of 
the clump formation. Jacquet & Krumholz (2011) performed 
linear stability analysis of a plane-parallel superposition of two 
media immersed in radiation field. They examined the stabil- 
ity in the optically thin and thick limit, finding that both cases 
could be Rayleigh-Taylor unstable. Such radiation-induced 
Rayleigh-Taylor instability also appears in massive star forma- 
tion system and H II region (see also Mathews & Blumenthal 
1977; Krumholz et al. 2009). 

The Rayleigh-Taylor instability, however, cannot explain the 
facts (2) and (3), which rather imply radiation processes being 
somehow involved. We thus further examine radiation-related 
instabilities. From Equation (14) and (15), radiation force 
could work as the positive feedback to grow the initial pertur- 
bation of the matter density, leading to a formation of inhomo- 
geneous density pattern. Shaviv (2001a) made a global linear 
stability analysis for the optically thick, radiation-dominated 
atmospheres of the super-Eddinton outflow from stars. He has 
clearly shown that perturbations with wavelengths on the order 
of pressure scale-height H grow on the order of sound cross- 
ing timescales under a fixed radiation-temperature condition at 
the bottom. This is a radiation-hydrodynamic instability and is 
expected to create inhomogeneous (porous) density structure 
just below photosphere of the objets (see also Shaviv 2001b; 
Fukue 2003). Note that one pressure scale-height roughly cor- 
responds to one optical depth just below the photosphere. That 
is, this instability may help forming clumps of one optical 
depth in the present study. 

We here found the anti-correlation between the matter den- 
sity and the absolute value of radiation force per unit mass in 
the present analysis. In contrast, Shaviv (2001a) claimed anti- 
correlation between matter density and radiation energy den- 
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R/r, 

Fig. 7. Two-dimensional structure of a supercritical accretion flow and 
the associated outflow simulated by the same radiation-MHD code but 
with no magnetic field. Green lines indicate the surface which upward 
radiation force equals downward gravity force of the central black hole. 
Similar clumpy outflow structure to that found in Figure 1 is confirmed. 

sity, which we do not confirm in our simulations. The reason 
for this should be clarified in a future work. 

In this simulation, we assumed that the flow is axisymmetric 
with respect to the rotation axis; i.e., the dense outflowing gases 
in the present simulation are not clumpy but annular. However, 
since the Rayleigh-Taylor and the radiation-hydrodynamic in- 
stability would work the clumpy structure formation even in 
the three-dimensional simulations. 

4.2. Role of Magnetic Field in Clump Formation 

When strong, ordered magnetic field are present in radiation- 
dominated atmosphere, the magnetic photon bubble instabil- 
ity may also work to produce inhomogeneous structure on the 
condition of Fq > E^Cs (Arons 1992; Gammie 1998; Blaes 
& Socrates 2003). Turner et al. (2005) performed the lo- 
cal radiation-MHD simulations of the standard thin (Shakura- 
Sunyaev) disk with the initial mass accretion rate 10% of the 
Eddington limit for a 10% radiative efficiency, and reported 
the occurrence of the magnetic photon bubble instability, giv- 
ing rise to inhomogeneous density structure. They also assert 
that the radiation-hydrodynamic instability of Shaviv (2001a) 
does not occur, although the this instability grows faster than 
the magnetic photon bubble instability, since the fixed radiation 
temperature lower boundary is not satisfied. 

To see if the existence of magnetic field is essential for 
the clumpy structure formation in our simulations, we per- 
formed non-magnetic radiation-hydrodynamic simulations. In 
this simulation, we adopted the data of the radiation-MHD 
calculation at the elapsed time of 7 s as the initial condition 
and then run the same code but by artificially vanishing mag- 
netic field everywhere. Once magnetic field are made zero, 
they never appear forever (see induction equation in Equation 
5). We should also note that there is no accretion motion in 
this simulation, since the disk viscosity, which is of a mag- 
netic origin, disappears. This simulation is nevertheless useful 



to see the dynamical properties of the non-magnetized outflow 
which takes place on much shorter timescales than the viscous 
timescale of the underlying disk. 

Figure 7 illustrates the resultant density distribution in the 
whole computational domain at the elapsed time of t = 9 s. 
The green lines indicate the surface of x^o/c = pIV-j/'pnI. It 
clearly shows again clumpy outflow structure above z ^ 250rs, 
where the upward radiation force overcomes the downward 
gravity force. We further confirm the typical clump size to 
be lOrs from the ACF analysis and anti-correlation between 
the matter density and the radiation force in the non-magnetic 
case, as well. We can thus conclude that magnetic field is not 
essential for creating clumpy outflows. In other words, the 
magnetic photon bubble instability is not a primary cause of 
clump formation. This result, however, does not exclude a 
possible occurrence of magnetic photon bubble instability in 
a real situation, since the wavelength with the fastest growth 
rate is ^ {^pgns/ Eo)H ^ H, which is too short to resolve by 
the present simulations. We need finer grid spacing, than the 
present case, to see the effects of magnetic photon bubble in- 
stability, if it really occurs. 

4.3. Comparison with the Observational facts 

As was mentioned in Introduction, the emergence of out- 
flow from various types of black hole objects has been reported 
through a number of X-ray observations. The clumpy nature 
of the outflow has also been suggested recently by the observa- 
tions of luminous accretion flows, such as ultra-luminous X-ray 
sources (ULXs) and bright AGNs. That is, the clumpy outflow 
nature does not depend on black hole masses observationally. 
We note that our results of the inflow-outflow structure also ap- 
ply to the cases with supermassive black holes, since Thomson 
scattering dominates over absorption opacities in the supercriti- 
cal accretion regimes, whatever black hole masses will be. The 
matter density decreases with an increase of the black hole 
mass, p oc Af for a fixed mass accretion rate. Hence, the 
radiation force on the Thomson scattering (oc p^) is dominant 
over the bound-free and free-free absorptions (oc p) even in the 
case of massive holes. 

We actually performed the same radiation-MHD simulations 
for the case with a supermassive black hole {M = 10^ Mq). 
The resultant density contours shown in Figure 8 clearly 
demonstrate the emergence of similar clumpy outflow patterns. 
By calculating ACFs, we also confirmed that the typical clump 
size is ^ lOrs, corresponding to Td = KosPci^ci ^ 1 (since 
Pel oc AI^^ while tc\ <x. M). 

Based on the variability study of the ULX, NGC 5408 X- 
1, Middle ton et al. (2011) found that it has a similar spectra 
with the black hole binary (BHB), GRS 1915H-105 and some 
extreme Narrow-Line Seyfert Is. They conclude that the un- 
derlying accretion flow should be supercritical and the associ- 
ated outflow has clumpy structure with the variability on time- 
scales of several tens of seconds (see also Middleton et al. 
2012). Tombesi et al. (2010) detected blueshifted Fe K absorp- 
tion lines from several AGNs whose the bolometric luminosi- 
ties are estimated to be close to the Eddington one, and sug- 
gested the presence of outflows from the central regions with 
mildly relativistic velocities, in the range 0.04-0. 15c: so called 
Ultra-Fast Outflows (UFOs). The ionization of the absorbers 
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Fig. 8. Two-dimensional structure of a supercritical accretion flow and 
the associated outflow around a black hole with M = 10*Mq at the 
elapsed time of t = 7 X lO'^ s. Green lines indicate the surface which 
upward radiation force equals downward gravity force of the central 
black hole. The clumpy structure does not depend on the black hole 
mass. 

is in the range ^ ~ 10^-10^ erg cm s~^, and the column den- 
sity is in the interval A^h 10^^-10^^ cm~^ (Tombesi et al. 
2011). The observed spectral variability on time-scales of ^ 
days (e.g., Braito et al. 2007; Cappi et al. 2009) suggests that 
the outflows have clumpy structure (Tombesi et al. 2012). 

Let us compare our results with these observational facts. It 
is important to note that once a clump passes across our line 
of sight, it will produce significant absorption, since the clump 
optical depth is about unity. To see how frequently such ob- 
scuration occurs, we calculate the spatial covering factor of the 
clumps, C, by integrating the fraction of sky covered by each 
clump element. The factor can be derived by 



ricX 



'sm0drded(p, (17) 



where rid is the number density of the clump, rph is the loca- 
tion of the photosphere, and {62 — di) is the opening angle of 
clumpy outflow, respectively (e.g., Bottorff & Ferland 2001). 
The spherical polar coordinates (r, 6, ip) are used. Note that 
the clump size £c\ is not radius but diameter The number den- 
sity is given by 



n-ci ' 



(18) 



where we use the following relation: the mass outflow rate 
Mout — pvr^; the matter density of the outflow p ^ mc\nc\', 
the unit mass of the clump md = 7r£^jPci/6; the solid angle of 
the outflow = 2 J^^ sm0d0dip. Using Equation (12), the 
covering factor is estimated to be 



C - 0.3t- 



IOLe/ 




Vr \ ( _[ph_\ 



(19) 



Note that the covering factor does not depend on the black hole 



mass nor the size of the clump, as long as the clump size is 
expressed in terms of the Schwarzschild radius. We thus expect 
occasional obscuration of the light from the center by clumpy 
outflow, being thus in agreement with the observations. 

Next, we estimate the photoionization parameter of clumps. 
The photoionization parameter, ^, is first proposed by Tarter et 
al. (1969) and is expressed by 

e = ^, (20) 

where Lx is the luminosity in the X-ray band. Using the rela- 
tion, Tci = nciCTT^cb the photoionization parameter is estimated 
to be 



^loV; 



Lx 
O.IXt 



VlOrs J VSOOrs 



erg cm s 



(21) 



Equation (21) indicates that clumpy outflows are mildly ion- 
ized (Kallman & McCray 1982). Although the parameter does 
not explicitly depend on the black hole mass, strictly speaking, 
the luminosity in the X-ray band Lx depends on the black hole 
mass; Lx ^ -^boi for Galactic sources and Lx ^ O.lLboi for 
AGNs. Here, Lboi is the bolometric luminosity of the source. 
We assume Lboi ~ -^e in Equation (21). We thus expect line 
absorption features to be observed, in agreement with the ob- 
servations. 

Finally, let us estimate the variability timescale of the 
clumpy outflow t^x, which is given by 



, "1,1 



4i 



(22) 



We here assume that the angular momentum, which was 
frph at the photosphere (at r = rph), is conserved even 
r. We thus estimate 



when it travels to the radius. 



tr. 



'6 



\W^MqJ \250rsJ 



500rs 



day, 



(23) 



in good agreement with the observations of AGN. This estima- 
tion cannot explain the long variation timescales of the ULXs, 
on the order of several tens of seconds, unless the black hole is 
large (Middleton et al. 2011). If the same mechanism applies 
to ULXs, they should contain supermassive black holes. 

The simulation results can account for the basic properties 
of the UFOs of AGN observed with X-ray. We should note, 
however, that we adopt gray approximation in the present sim- 
ulations, meaning that we cannot properly discuss the spectral 
energy distribution (SED). More detailed comparison with the 
observations, in terms of absorption spectra, is left for future 
issue. 

4.4. The Origin of Broad-Line Clouds? 

The broad-line region (BLR) is a standard AGN ingredient 
(e.g. Ferland et al. 1992; Peterson 1997): it is located inside the 
obscuring torus, above and below the accretion flow, and at a 
distance of parsec scale from the supermassive black hole. Can 
we explain BLR clouds by clumpy outflow? From simulation 
data, the electron number density and the gas temperature of 
the clumps are estimated ~ Pci/'m-p ~ 10^^ cm^'^ and ^ lO'* 



10 



S. Takeuchi, K. Ohsuga, and S. Mineshige 



[Vol. , 



K for M = 10''A/q, since the gas temperature is proportional 
to M^'^/^ for radiation-dominated atomosphere. The volume 
filling factor is expressed by 

^ iVei«i/6) _ N^y^l 

where A'ci = /// ricir^ sin 9drd9d(p is the number of clumps, 
and rsLR is the size of the BLR. By using some typical values, 
we got the filling factor 

/ Mout \ ( ""r \ ( ^cl \ ( ^BLR \ 

yioLE/cy [oTcJ [To^j l^ioVs J ■ ^^^^ 

This value is consistent with the observational values of the 
BLR clouds (Peterson 1997). In this estimation, we assumed 
that each clump keeps its form up to the distance of tblr = 
lO^rs ~ Ipc. We check if clumps survive by comparing var- 
ious timescales. The gas temperature of the clumps hardly 
change due to the matter-radiation coupling occur in them, 
whereas that of the ambient media would remain to be hot since 
the radiative cooling cannot be effective. The thermal conduc- 
tion time scale, which is expressed by SugpUcikBi'^iT"^^^ , is 
also much lager than the dynamical time scale, where Ksp is 
the Spitzer conductivity coefficient. Thus the clumps are ex- 
pected to remain to large domain. To examine the stability 
of the clumps, the study with more larger simulation box is 
needed. 

The estimated typical values from the super- Eddington out- 
flow with the mass outflow rate of A/out = IOLe/c^ is in agree- 
ment with the observational facts. Importantly, the photon lu- 
minosity of the underlying accretion flow is ~ in spite of the 
large mass accretion and outflow rate, since the a large amount 
of photons inside the accretion flow is trapped and swallowd 
to the central black hole (so-called the photon-tapping effect) 
(Ohsuga et al. 2005). 

Elitzur (2012) attempts to explain the observations of type I 
and II AGNs based on the assumption that clumpy outflow gas 
originating from the accretion flow is distributed around lumi- 
nous AGN. His model can nicely explain the observed BLR 
disappearance at low luminosity (Elitzur & Ho 2009). 

4.5. Possible line-driven outflow 

Although we only consider the continuum radiation force 
in the present study, line force may also work in luminous 
black hole inflow-outflow system. The line-driven force has 
been suggested to the one of the most promising mechanisms 
of AGN outflow, since this can explain both acceleration and 
ionization states of the outflow materials (Proga et al. 2000). 
When matter absorbs UV photons, it gets momentum from 
them, and is accelerated towards the direction opposite to that 
the photons come from. The ionized X-ray photon distribution 
determines the outflow structure (Castor et al. 1975; Risaliti 
& Elvis 2010; Nomura et al. 2012). Since the SED of the 
underlying accretion flow strongly depends on the dynamics 
of line-driven outflows, it is unclear whether the line-driven 
mechanism does really work or not in the present study (note 
that we assume the gray approximation for radiation transfer). 
Kurosawa & Proga (2009) performed the three-dimensional 



hydrodynamical simulations of line-driven wind, reporting the 
clumpy outflow structure due to the rotation shear. Their find- 
ing is, however, independent of our work, in which radiation- 
hydronamical effects seem to play a critical role. The per- 
formance of the global radiation-MHD simulation of the line- 
driven outflow is future work. 

5. Conclusion 

In this paper, we performed the two-dimensional global 
radiation-MHD simulations of supercritical accretion flows 
onto black holes in the larger computational box than our pre- 
vious one and examined the properties of the outflow from the 
accretion flow in details. Here are our new findings: 

• The outflows associated with supercritical (or super- 
Eddington) accretion flows have a clumpy structure 
above heights of ~ 2507's. The typical clump size is 
^ lOrs, which corresponding to about one optical depth, 
and their shapes are slightly elongated along the outflow 
direction. In the clumpy outflow region, a clear anti- 
correlation is seen between matter density and the ab- 
solute value of radiation force per unit mass. 

• Rayleigh-Taylor instability is the most plausible cause 
of the clump formation, since the clumpy structure ap- 
pears in the layer where upward radiation force over- 
comes downward gravity force. In addition, a radiation- 
hydrodynamic instability, which arises when radia- 
tion funnels through radiation-pressure supported atmo- 
sphere, may also help forming clumps of one optical 
depth. Magnetic filed is not essential for this clumpy 
structure formation. 

• The spatial covering factor of the clumps are estimated 

0.3 for typical parameters of supercritical accretion 
flow, regardless of the black hole mass nor the size of 
the clump. This mean that presence of the occasional 
obscuration of the light from the center by clumpy out- 
flow. For the clumpy outflows from supermassive black 
hole with Al = IO^A/q, the photoionization parameter 
and the variability timescale of the clump are estimated 
~ 10^ erg cm s^^ and ^ 6 day, in good agreement with 
the observational results. 

• If the clumps remain to parsec scale, the volume filling 
factor is estimated ~ 10^^, and the clumpy outflow are 
consistent with the observational values of BLR clouds. 
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Appendix. Correlation Function 

The correlation function is a measure of similarity of two 
waveforms as a function of a lag. Correlation analysis is used 
to find periodic patterns and/or coherent lengths in data. The 
correlation function /con (i) of two sample populations x and 
y as a function of the lag L is calculated by, 



{xk+L-x){yk-y) 



k=0 



/cori(i) — < 




for L< 



*;=0 



N-1 



fc=0 



■Af-1 



for L> 



fc=0 



where x and y are the means of the sample populations x = 
{xo,xi,X2,..;Xn-i) and y ^ (?/o, 2/1, J/2, J/w-i), respec- 
tively. The correlation function is called auto-correlation func- 
tion (ACF) in the case of same waveform (x = y), and cross- 
correlation function (CCF) in the case of different waveforms 
(x + y)- 
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